Exploración de las relaciones entre aguas subterráneas y variables del suelo

Contexto

Los humedales altoandinos son piezas clave de los balances de agua, energía y carbono en la región del Altiplano, procesos que dependen de la disponibilidad de aguas subterráneas cercanas a la superficie. Pese a su relevancia, estos ecosistemas enfrentan amenazas crecientes debido al cambio climático y la actividad minera. Además, los métodos de monitoreo hidrogeológico actuales son caros y difíciles de implementar en estos ambientes remotos.

El objetivo de este trabajo es explorar las relaciones entre profundidad del agua subterránea y variables medidas en el suelo para evaluar el desarrollo de metodologías de monitoreo alternativas, usando el Salar del Huasco como caso de estudio.

Figura 1. Localización del Salar del Huasco y sus estaciones de monitoreo

En la actualidad, se cuenta con dos estaciones de monitoreo en el sitio: norte y sur. En este trabajo se emplearán los datos recolectados en la estación norte, donde se cuenta con un piezómetro y tres sensores de suelo ubicados a diferentes profundidades (15, 30 y 48 cm). El piezómetro reporta la profundidad del nivel freático (m), mientras que los sensores de suelo, datos de contenido de agua (m³/m³)y temperatura (°C).

Análisis exploratorio

Los datos crudos entregados por los instrumentos fueron preprocesados previo a este análisis. Los scripts correspondientes están disponibles aquí y el d. Los datos crudos y procesados están disponibles aquí. El código completo asociado a las visualizaciones presentadas está disponible en este notebook.

Series temporales

En primer lugar, se graficaron las series de tiempo de profundidad del nivel freático, así como de contenido de agua y temperatura del suelo a las tres profundidades disponibles (Figura 2). Cabe mencionar que el eje Y del primer panel está invertido para facilitar la interpretación de la profundidad de las aguas subterráneas.

Se observa que la profundidad del nivel freático es variable en el tiempo, alcanzando máximos en primavera y mínimos en invierno. Aparentemente, el contenido de agua en el suelo a 15 y 30 cm sigue un patrón similar al de las aguas subterráneas, puesto que comparten momentos de ascenso y descenso. El contenido de agua a 48 cm no muestra mayor variabilidad durante el año. Por otro lado, la temperatura del suelo tiene un marcado ciclo estacional, con valores máximos en verano y mínimos en invierno, por lo que, no se ve relacionada con la profundidad del nivel freático a primera vista.

# Creacion de la figura base
time_series = sns.relplot(
    data=long_df,
    x='Timestamps',
    y='Value',
    hue ='Depth-label',             # Diferencia colores segun profundidad
    row='Variable-label',           # Una fila por cada variable
    palette=palette_dict,
    hue_order=full_depth_order,
    kind='line',                    # Grafico de lineas
    height=3,
    aspect=2.5,
    facet_kws={'sharey': False, 'sharex':True}      # Eje Y independiente, eje X compartido
)

# Eliminacion de elementos incluidos por defecto
time_series.set_titles('')
time_series.set_xlabels('')
time_series._legend.remove()

# Iteracion sobre facets para aplicar estilos particulares
for i, (ax, title) in enumerate(zip(time_series.axes.flat, time_series.row_names)):

    # Etiquetado de los ejes Y de cada facet
    ax.set_ylabel(title, fontweight='normal')

    # Formateo de fechas en el eje X
    date_format = mdates.DateFormatter('%b %Y')
    ax.xaxis.set_major_formatter(date_format)

    # Configuracion particular para el primer facet
    if 'Prof. nivel freático (m)' in title:
        
        # Invierte el eje Y
        ax.invert_yaxis()

        # Traza profundidad de sensores de suelo como lineas horizontales
        sensor_depths_m = [0.15, 0.30, 0.48]
        for depth, label in zip(sensor_depths_m, soil_depth_order):
            ax.axhline(y=depth, color=palette_dict[label], linestyle='--', linewidth= 0.9, alpha=0.8)

# Eliminacion de la leyenda de objetos asociados al piezometro
handles, labels = time_series.axes[0,0].get_legend_handles_labels()
clean_handles = [h for h, l in zip(handles, labels) if l != 'NA']
clean_labels = [l for l in labels if l != 'NA']

# Creacion de nueva leyenda
time_series.figure.legend(
    handles=clean_handles,
    labels=clean_labels,
    loc='lower center',
    bbox_to_anchor=(0.5, 0),
    ncol=3,
    title="Profundidad de los sensores de suelo",
    frameon=False
)

# Establecimiento de titulo
# time_series.figure.suptitle("Evolución de variables piezométricas y edáficas",
#                             x=0.5, y=0.98, fontsize=14)

# Ajustes de espaciado
plt.tight_layout()
plt.subplots_adjust(top=0.93, bottom=0.11) # Ajusta espacio disponible para facets

plt.show()

Figura 2. Series temporales de profundidad del nivel freático, contenido de agua y temperatura del suelo.

Correlaciones cruzadas

Para evaluar la existencia de posibles desfases temporales entre las series de datos, se ejecutó un análisis de correlaciones cruzadas. En él, se calculó la correlación de Pearson entre la serie de profundidad del nivel freático y las series de contenido de agua y temperatura del suelo, desplazadas temporalmente hasta ±15 días (Figura 3). Se definió esta temporalidad porque no son esperables relaciones significativas en periodos mayores. Previo al análisis, las series fueron procesadas para alcanzar la estacionariedad.

Para la variable de contenido de agua, se observa que la mayor correlación con el nivel freático se da en el desfase cero, y es negativa para los tres sensores. Esto indica que un aumento del contenido de agua se acompaña de una disminución de la profundidad de las aguas subterráneas, y viceversa. La relación más fuerte se da en el sensor a 30 cm. Respecto a la temperatura del suelo, se constatan relaciones inversas en el desfase -1 en los sensores a 30 y 48 cm, lo que significa que un aumento en la temperatura suele anteceder en un día un aumento en la profundidad del nivel freático en los sensores más profundos.

# Creacion de la figura base
cross_corr = sns.relplot(
    data=cross_corr_df,
    x='Lag',
    y='Cross-correlation',
    hue='Depth-label',
    row='Variable-label',
    palette=palette_dict,
    kind='line',
    height=3,
    aspect=2.5
)

# Eliminacion de elementos incluidos por defecto
cross_corr.set_titles('')
cross_corr.set_xlabels('')
cross_corr._legend.remove()

# Iteracion sobre facets para etiquetar los ejes Y
for ax, title in zip(cross_corr.axes.flat, cross_corr.row_names):
    ax.set_ylabel(title)

# Etiquetado del eje X
ax.set_xlabel('Días de desfase')
# Trazado de lineas de referencia
for ax in cross_corr.axes.flat:

    # Trazado de lineas 0, 0
    ax.axvline(0, color='gray', linestyle='--', linewidth=0.8, alpha=0.6)
    ax.axhline(0, color='black', linestyle='-', linewidth=1, alpha=0.7)

    # Trazado de limites de confianza
    ax.axhline(conf_cross_corr, color='gray', linestyle='--', linewidth=0.8, alpha=0.5)
    ax.axhline(-conf_cross_corr, color='gray', linestyle='--', linewidth=0.8, alpha=0.5)
    ax.axhspan(-conf_cross_corr, conf_cross_corr, color='gray', alpha=0.05)

    # Anotacion del limite de confianza
    ax.text(
        x=0.82,
        y=0.87,
        s='Intervalo de confianza al 95%',
        transform=ax.transAxes,
        fontsize=8,
        color='gray'
    )  

# Creacion de nueva leyenda
cross_corr.figure.legend(
    loc='lower center',
    bbox_to_anchor=(0.53, -0.07),
    ncol=3,
    title="Profundidad de los sensores de suelo",
    frameon=False
)

# Establecimiento de titulo
# cross_corr.figure.suptitle("Correlaciones cruzadas entre nivel freático y variables del suelo", 
#                            x=0.5, y=0.98, fontsize=14)

# Ajustes de espaciado
plt.tight_layout()
plt.subplots_adjust(top=0.93, bottom=0.11) # Ajusta espacio disponible para facets

plt.show()

Figura 3. Correlaciones entre series de profundidad del nivel freático y variables del suelo, desfasadas temporalmente.

Gráficos de dispersión

A fin de identificar las relaciones entre nivel freático y variables del suelo, se generaron gráficos de dispersión a distintas profundidades con datos diferenciados por mes, para visualizar su evolución temporal (Figura 4). Notar que los ejes Y están invertidos en todos los gráficos.

En términos generales, en los sensores a 15 y 30 cm se observa una relación inversa entre la profundidad del nivel freático y el contenido de agua del suelo, que es aparentemente no lineal. Por otra parte, en el sensor a 48 cm no existe una variación importante del contenido de agua, independientemente de la profundidad del nivel. Del mismo modo, la profundidad del agua subterránea parece ser independiente de la temperatura, siendo las variaciones atribuibles a los cambios de estación.

# Creacion de la figura base
scatter_plot = sns.relplot(
    data=scatter_df,
    x='Value',                          # Eje X: variables de suelo
    y='Groundwater-depth-m',            # Eje Y: profundidad nivel freatico
    hue='Month',                        # Coloreado por mes
    row='Depth-label',                  # Una fila por profundidad
    col='Variable-label',               # Una columna por variable
    row_order=soil_depth_order,
    palette='twilight_shifted',         # Paleta ciclica
    kind='scatter',
    height=3,
    aspect=1,
    facet_kws={'sharey': True, 'sharex': 'col'},    # Eje Y compartido, Eje X: compartido por columnas
    s=50,
    alpha=0.9,
    legend='full'
)

# Eliminacion de elementos incluidos por defecto
scatter_plot.set_titles('')
scatter_plot.set_xlabels('')
scatter_plot.set_ylabels('')

# Etiquetado del eje X
scatter_plot.axes[-1, 0].set_xlabel("Contenido de agua (m³/m³)")
scatter_plot.axes[-1, 1].set_xlabel("Temperatura (°C)")

# Inversion del eje Y (+ margen)
y_min, y_max = scatter_df['Groundwater-depth-m'].min(), scatter_df['Groundwater-depth-m'].max()
padding = (y_max - y_min) * 0.05
scatter_plot.set(ylim=(y_max + padding, y_min - padding))

# Etiquetado del eje Y
scatter_plot.figure.supylabel(
    'Profundidad del nivel freático (m)',
    x=0.04,
    size=plt.rcParams['axes.labelsize'], # Igualar tamaño al eje X
    )

# Titulo de leyenda
legend = scatter_plot._legend
legend.set_title('Mes')

# Etiquetado de meses en leyenda
month_names = {str(i): calendar.month_abbr[i] for i in range(1, 13)}
for t in legend.texts:
    if t.get_text() in month_names:
        t.set_text(month_names[t.get_text()])

# Iteracion sobre los facets para agregar etiqueta de profundidad
for ax, label in zip(scatter_plot.axes[:, -1], soil_depth_order):
    ax.text(
        x=0.95,
        y=0.95,
        s=f'{label}',
        transform=ax.transAxes,
        ha='right',
        va='top',
        fontsize=10,
    )

# Establecimiento de titulo
# scatter_plot.figure.suptitle("Relación entre variables de suelo y profundidad del nivel freático",
#                              x=0.5, y=0.98, fontsize=14)

# Ajustes de espaciado
plt.tight_layout()
plt.subplots_adjust(top=0.94, right=0.85) # Ajusta espacio disponible para facets

plt.show()

Figura 4. Relaciones entre profundidad del nivel freático y variables del suelo diferenciadas mensualmente. Eje Y invertido

Correlaciones móviles

Con el objetivo de visualizar la evolución de las correlaciones entre variables en diferentes estaciones del año, se hizo un análisis de correlación empleando ventanas móviles de 61 días. En él, se calculó la correlación de Pearson entre las series estacionarias para cada día disponible usando las observaciones de ±30 días, además del día presente (Figura 5). Se definió una ventana de dos meses para contar con un número de observaciones suficiente para el análisis estadístico, además de suavizar los resultados obtenidos y facilitar la interpretación.

Las relaciones más fuertes entre contenido de agua y profundidad del nivel freatico se dan durante el verano y son inversas, siendo las mayores a la profundidad de 30 cm. Durante la primavera se da una relación significativa directa en los sensores a 15 y 30 cm, lo que indica que un aumento del contenido de agua se acompaña de un aumento de la profundidad del nivel. En la variable temperatura se da un patrón opuesto: en primavera se dan relaciones inversas, mientras que en verano son directas.

# Creacion de la figura base
roll_corr = sns.relplot(
    data=rolling_corr_df,
    x='Timestamps',
    y='Rolling-correlation',
    hue='Depth-label',
    row='Variable-label',
    palette=palette_dict,
    kind='line',
    height=3,
    aspect=2.5
)

# Eliminacion de elementos incluidos por defecto
roll_corr.set_titles('')
roll_corr.set_xlabels('')
roll_corr._legend.remove()

# Iteracion sobre facets para etiquetar los ejes Y
for ax, title in zip(roll_corr.axes.flat, roll_corr.row_names):
    ax.set_ylabel(title)

# Formateo de fechas en el eje X
    date_format = mdates.DateFormatter('%b %Y')
    ax.xaxis.set_major_formatter(date_format)

# Trazado de lineas de referencia
for ax in roll_corr.axes.flat:

    # Trazado de linea 0
    ax.axhline(0, color='black', linestyle='-', linewidth=1, alpha=0.5)

    # Trazado de limites de confianza
    ax.axhline(conf_rolling_corr, color='gray', linestyle='--', linewidth=0.8, alpha=0.5)
    ax.axhline(-conf_rolling_corr, color='gray', linestyle='--', linewidth=0.8, alpha=0.5)
    ax.axhspan(-conf_rolling_corr, conf_rolling_corr, color='gray', alpha=0.1)

    # Anotacion del limite de confianza
    ax.text(
        x=0.82,
        y=0.8,
        s='Intervalo de confianza al 95%',
        transform=ax.transAxes,
        fontsize=8,
        color='gray'
    )    

# Creacion de nueva leyenda
roll_corr.figure.legend(
    loc='lower center',
    bbox_to_anchor=(0.5, -0.03),
    ncol=3,
    title="Profundidad de los sensores de suelo",
    frameon=False
)

# Establecimiento de titulo
# roll_corr.figure.suptitle(f'Correlación entre nivel freático y variables del suelo en ventanas móviles de {window_days} días', 
#                            x=0.5, y=0.98, fontsize=14)

# Ajustes de espaciado
plt.tight_layout()
plt.subplots_adjust(top=0.93, bottom=0.11) # Ajusta espacio disponible para facets

plt.show()

Figura 5. Evolución temporal de las correlaciones entre profundidad del nivel freático y variables del suelo en ventanas móviles de dos meses.

Diagrama de Hövmoller

Finalmente, se construyó un diagrama de Hövmoller para representar la evolución temporal de las variables de suelo en un perfil vertical, al cual se superpuso la profundidad de las aguas subterráneas (Figura 6). El perfil se obtuvo interpolando cada un centímetro los datos entre los sensores de suelo, para cada día disponible.

En el diagrama superior se ve cómo la profundidad de las aguas subterráneas tiene un efecto sobre la distribución vertical del agua en el suelo. Por el contrario, el diagrama de temperaturas no muestra una relación aparente con el nivel freático, sin embargo, permite visualizar que la temperatura en el perfil de suelo está condicionada a lo que ocurre en superficie. Durante el invierno es el enfriamiento superior el que se transmite hacia capas profundas, mientras que en verano la existencia de pulsos de calor contribuyen a aumentar la temperatura en profundidad. La distribución temporal de los pulsos sugiere que se trata de eventos de precipitación.

# Definicion de funcion para generar un mapa de calor
def heatmap(data, **kwargs):
    
    # Pivoteo de los datos largos a matriz. La funcion de mapeo de heatmap
    # de matplotlib (pcolormesh) requiere datos en formato ancho
    pivot_data = data.pivot(index='Depth_m', columns='Date', values='Value')

    # Defincion de objetos que guardan los elementos de la matriz
    X = pivot_data.columns  # Fechas
    Y = pivot_data.index    # Profundidades
    Z = pivot_data.values   # Valores

    # Asignacion condicional de colores y etiquetas segun variable
    variable = data['Variable-label'].iloc[0]
    if 'Temperatura (°C)' in variable:
        palette = 'OrRd'
        label = 'Temperatura (°C)'
    else:
        palette = 'PuBu'
        label = 'Contenido de agua (m³/m³)'

    # Renderizado de la grilla interpolada usando pcolormesh
    mesh = plt.pcolormesh(X, Y, Z, cmap=palette, shading='auto')

    # Definicion del objeto ax
    ax = plt.gca()

    # Trazado de nivel freatico en los axes
    gw_df = kwargs.get('groundwater_data')
    col_piezo = 'Piezometer_NA_groundwater-depth_m' 
    ax.plot(
        gw_df.index, 
        gw_df[col_piezo], 
        color='black',
       linestyle='-',
        linewidth=1.5,
        label='Nivel freático'
        )

    # Creacion de leyenda
    cbar = plt.colorbar(mesh, 
                 ax=ax, 
                 label=label, 
                 pad=0.02, 
                 aspect=10,
                 shrink=0.7)
    cbar.set_label(
        label, 
        rotation=-90,
        labelpad=15
        )

hovmoller = sns.FacetGrid(
    data=hovmoller_df,
    row='Variable-label',
    height=4,
    aspect=2.2,
    sharex=True,
    sharey=True
)

# Mapeo del heatmap
hovmoller.map_dataframe(heatmap, groundwater_data=wide_df)

# Eliminacion de elementos incluidos por defecto
hovmoller.set_titles('')

# Iteracion sobre facets
for i, ax in enumerate(hovmoller.axes.flat):
    
    # Añadido de leyenda del nivel freatico 
    if i == 1:
        ax.legend(loc='upper right', frameon=False, fontsize='small', framealpha=0.8)
    
    # Formateo de fecha
    date_fmt = mdates.DateFormatter('%b %Y')
    ax.xaxis.set_major_formatter(date_fmt)

    # Añadido de profundidad de sensores
    for d in [0.15, 0.30, 0.48]:
        ax.axhline(d, color='gray', linestyle='--', alpha=0.3, linewidth=1)
        ax.text(wide_df.index[0], d, f' {d*100:.0f} cm', color='gray', va='bottom', fontsize=8, alpha=0.7)

# Inversion del eje Y
y_min, y_max = 0, 0.7
hovmoller.set(ylim=(y_max, y_min))

# Etiquetado del eje Y
hovmoller.figure.supylabel(
    'Profundidad del perfil (m)',
    x=0.04,
    size=plt.rcParams['axes.labelsize'], # Igualar tamaño al eje X
    )

# Establecimiento del titulo
# hovmoller.figure.suptitle("Variables de suelo interpoladas en un perfil vertical",
#                             x=0.5, y=0.98, fontsize=14)

# Ajustes de espaciado
plt.tight_layout()
plt.subplots_adjust(top=0.94, bottom=0.11, left=0.1) # Ajusta espacio disponible para facets

plt.show()

Figura 6. Dinámica temporal de las variables de suelo en el perfil con datos interpolados verticalmente.

Desafíos

Desarrollar una conclusión a partir de los resultados iniciales y cuáles serían próximos pasos a desarrollar para continuar con la investigación.